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Abstract 

The relation between the results of cosmological N-body simulations, and the continuum theoretical 
models they simulate, is currently not understood in a way which allows a quantification of N 
dependent effects. In this first of a series of papers on this issue, we consider the quantification of 
such effects in the initial conditions of such simulations. A general formalism developed in jlj allows 
us to write down an exact expression for the power spectrum of the point distributions generated 
by the standard algorithm for generating such initial conditions. Expanded perturbatively in the 
amplitude of the input (i.e. theoretical, continuum) power spectrum, we obtain at linear order 
the input power spectrum, plus two terms which arise from discreteness and contribute at large 
wavenumbers. For cosmological type power spectra, one obtains as expected, the input spectrum 
for wavenumbers k smaller than that characteristic of the discreteness. The comparison of real space 
correlation properties is more subtle because the discreteness corrections are not as strongly localised 
in real space. For cosmological type spectra the theoretical mass variance in spheres and two point 
correlation function are well approximated above a finite distance. For typical initial amplitudes this 
distance is a few times the inter-particle distance, but it diverges as this amplitude (or, equivalently, 
the initial red-shift of the cosmological simulation) goes to zero, at fixed particle density. We discuss 
briefly the physical significance of these discreteness terms in the initial conditions, in particular 
with respect to the definition of the continuum limit of N-body simulations. 

PACS numbers: 98.80.-k, 05.70.-a, 02.50.-r, 05.40.-a 



I. INTRODUCTION 

The goal of dissipationless cosmological N-body simu- 
lations (NBS) is to trace the evolution of the clustering of 
matter under its self-gravity, starting from a cosmologi- 
cal time at which perturbations to homogeneity are small 
at the physically relevant scales (for reviews, see e.g. 
H [|, Qj). A fundamental question about such simulations 
concerns the degree to which they reproduce the evolu- 
tion of the simulated models. The problem arises because 
the N-body technique, which simulates a large number 
of particles evolving under their self-gravity (with some 
small scale regularization), is not a direct discretization 
of the theoretical models. The N-body approach is taken 
because it is not numerically feasible to simulate, at a 
useful level of resolution, the continuum Vlasov-Poisson 
equations describing the evolution theoretically. 

Since the Vlasov-Poisson system corresponds to an ap- 
propriately defined N — > oo of this particle dynamics 
H ||, NBS can be considered as related to the mod- 
els in this limit. The problem of discreteness is thus that 
of the relation of the results obtained from these sim- 
ulations, for typical statistical quantities characterising 
clustering, with those which would be obtained with such 
a simulation done with N — > oo particles. The existing 
studies of this "convergence" problem in the literature 
(e.g. |^|, [?], |8|, [I |ll)) are almost exclusively numerical, 



and consider the stability of different measured quanti- 
ties as a function of N. In the absence, however, of any 
analytic understanding of the possible N dependence of 
the results, such studies, which extend over a very mod- 
est range of N, cannot be conclusive. Different groups 
of authors have in fact drawn very different conclusions 
about the correctness of results for standard quantities 
at smaller scales. Some Q [l^] even place in question the 
validity of results for clustering amplitudes below the ini- 
tial interparticle distance, while such results are widely 
interpreted as physical in almost all current simulations. 
Further it is not specified in such studies how precisely 
the limit of large N should be taken, i.e., which other pa- 
rameters (e.g. box-size, force softening, initial red-shift) 
should be kept fixed or varied. These questions are be- 
coming of ever greater practical importance as the quan- 
tification of the precision of results from simulations is 
essential in order to confront cosmological models with a 
rich host of observations (see, e.g., fll3|). 

In this paper we address only one simple aspect of this 
problem: the relation between the discretized initial con- 
ditions (IC) of an NBS and the IC of the corresponding 
theoretical model. More specifically we study and quan- 
tify analytically the differences between the two for the 
two point correlation properties, in real space and re- 
ciprocal space, in the infinite volume limit (at a fixed 
particle density). The discreteness effects, i.e., the dif- 
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ferences between the continuous theoretical IC and the 
discrete IC of the actual NBS are then terms which de- 
pend on the particle density. We study these terms and 
their relative importance for different theoretical IC. We 
also study and check our analytic results with the aid of 
numerical simulations, performed in one dimension be- 
cause of the modest numerical cost of calculating the en- 
semble average using a large number of realisations. We 
underline that our analytic results apply to the infinite 
volume limit, i.e., to the case that the size of the simula- 
tion box goes to infinity at fixed particle density. They 
thus do not include effects associated with the finite size 
of the box. Such effects, which are quite distinct from 
those studied here, have been studied extensively else- 
where (see e.g. @, [H| H, |r| [l| 0), both analytically 
and numerically. 

The direct motivation for this work on IC in NBS came 
originally from some numerical studies of these issues by 
two groups |^0[ ^l], |22|, |22j ]. who have drawn quite dif- 
ferent conclusions about the accuracy with which the IC 
produced by the canonically used algorithm in NBS rep- 
resent the input IC . The results in this paper, which are 
essentially analytic, clarify the issues underlying the dis- 
cussion in and differences between these numerical stud- 
ies. Our conclusions are consistent with the findings of 
both sets of authors, and explain the differences. In short 
the authors of 0, are correct that certain real space 
properties, notably the mass variance in spheres, are in 
fact reasonably well represented for typical IC in NBS. 
The authors of pC] , however, are correct in diagnosing 
the important systematic differences between the actual 
and theoretical correlation properties in real space. In- 
deed one of our main findings is that there is a very non- 
trivial difference between the two spaces: while the dis- 
creteness of the underlying particle distribution is strictly 
localized in reciprocal space, this is not the case in real 
space. The result is that, in the limit of low amplitude 
initial density fluctuations — or, equivalently, high initial 
redshift for the simulation — the correlation properties of 
the input theoretical model are approximated well only 
in reciprocal space. Taking instead the limit that the 
particle density goes to infinity at fixed amplitude, the 
theoretical correlation properties are recovered in both 
real and reciprocal space, provided an appropriate cut- 
off is imposed at large wavenumber in the input PS. 

The implications of these results for what concerns 
the agreement between an evolved NBS and the evolved 
theoretical model of which it is the discretization is be- 
yond the scope of this paper. In a forthcoming paper 
p4j| , treating discreteness effects up to shell crossing, we 
will see that the evolution of an NBS deviates arbitrar- 
ily from its continuum counterpart as the initial red- 
shift increases at fixed particle density, while keeping the 
amplitude fixed at increasing particle one approximates 



increasingly well the continuum (fluid limit) evolution. 
Thus the results we find for the initial conditions here 
do indeed turn out to have physical significance for the 
question of discreteness effects in the evolved simulations. 

The algorithm used to generate IC in NBS which we 
analyze is in fact well defined, in the infinite volume limit, 
only for a certain range of asymptotic behaviors of the in- 
put theoretical PS. We specify here carefully these limits. 
While it turns out that they are not particularly relevant 
to current cosmological models, they are of interest in 
other contexts in which this algorithm may be used, no- 
tably in the study of gravitational clustering from other 
classes of initial conditions (e.g. ^||, ^6|). These prop- 
erties are also of interest in the context of statistical 
physics, where the problem of "realizability" of point pro- 
cesses is studied (see e.g. [^7], ^9|). Specifically we 
find that the algorithm has interesting limits for the case 
of very "blue" input PS: for the case of spectra with a 
small k behavior proportional to fc™, and n > 1, the real 
space variance is never that of the input model, while for 
n > 4 the reciprocal space representation is never faithful 
either. 

The paper is organized as follows. In the section which 
follows we briefly review the standard method for setting 
up IC for cosmological simulations, using the Zeldovich 
approximation. This also sets conventions for notation in 
the rest of the paper. In Sect. Ill we analyze the PS of the 



configurations of points generated in this way, comparing 
it with the PS of the input theoretical model. To obtain 
these results we use a very general exact result derived in 
0, which gives the two point properties of a point dis- 
tribution generated by superimposing an arbitrary corre- 
lated displacement field on an arbitrary initial stochastic 
point distribution. In the following section we consider 
how these properties described in /c-space translate into 
real space. Specifically we present a general qualitative 
analysis of the relation between the two point correlation 
properties of the IC and those of the input models. We 
treat specifically the mass variance in spheres, and the re- 
duced two point correlation function. For the latter case 
the comparison of the theoretical model and full IC is 
more difficult, because of the complicated non-monotonic 
form of the correlation function of the underlying point 
distributions. In the following section we illustrate, and 
verify our results quantitatively, using one dimensional 
numerical simulations. We choose to work in one dimen- 
sion for numerical economy, and because all the perti- 
nent questions can be posed equally well and answered 
in this case 2 . In the final section we summarize our re- 
sults, discuss what conclusions can be drawn concerning 
the papers mentioned above which motivated the present 
study, and finally briefly comment on the physical signif- 
icance of our results, which will be further developed in 



1 A more detailed account is given in the conclusions section below. 



2 The same is evidently not true when considering discreteness 
effect in the dynamics. 
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the companion paper [B4| . Several technical details in the 
paper, notably concerning the perturbative expansion of 
the exact expression for the PS of the generated IC, are 
given in three appendices. 



II. GENERATION OF IC USING THE 
ZELDOVICH APPROXIMATION 

The method which is used canonically for the genera- 
tion of IC in cosmological NBS is based on the so-called 
Zeldovich approximation (ZA) f30|| . It may be derived 
at linear order in a perturbative treatment of the equa- 
tions describing a self-gravitating fluid in the Lagrangian 
formulation pl| . It relates the initial position q of a fluid 
element to its final position 3 x through the expression 



x(q,t) =q+/(t)u(q), 



(1) 



i.e., it expresses the displacement of a particle as a sepa- 
rable function of the initial position q and the time t. The 
function f(t) is equal, up to an arbitrary normalization, 
to the growth factor of density perturbations derived in 
linear perturbation theory (see below). The vector field 
u(q) is thus proportional to both the velocity and accel- 
eration of the fluid element, and with a suitable normal- 
ization it can thus be taken to satisfy 



u(q) = -V • $(q) 



(2) 



where ^(q) is the gravitational potential at the initial 
time created by the density fluctuations. To set up IC 
representing a density field, one thus simply determines 
the associated potential through the Poisson equation 
and infers the appropriate displacements [and velocities, 
given f(t)] using Eq. (|l|) to apply to a set of points rep- 
resenting the unperturbed fluid elements. 

We note that to understand this algorithm for gener- 
ating IC representing a continuous density field it is not 
in fact necessary to invoke the ZA, nor anything which 
has specifically to do with gravity. These latter are rele- 
vant only for the determination of the velocity field. The 
only relation needed is in fact the continuity equation, 
which relates the velocity (and thus displacement) field 
in a fluid to the density fluctuations. At leading order in 
the density fluctuations it gives 



Sp(x) = -V • u(x). 
where the density fluctuation Sp(x) is defined by 



(3) 



(4) 



p(x) is the (continuous) density field, po the average den- 
sity, and u(x) is the displacement field. By inversion one 
can determine a displacement field which gives a desired 
density field. If one assumes further that the former is 
curl-free, and thus derivable as the gradient of a scalar 
field, one obtains a unique prescription for the displace- 
ment field which is identical to that given by the ZA as 
described above. 

In cosmological models the starting point for IC is not 
a specific density field, but a power spectrum (PS) of 
density fluctuations. The latter is defined as 



P(k) = lim 

V—>cx 



(\smi 

V 



(•5) 



where (...) denotes the average over an ensemble of re- 
alizations and Sp{k) denotes the Fourier transform (FT) 
of p(x) defined as 



6p(k) = / d d x6p(x)e- lU x . 
Jv 

It follows then from Eq. (||) that 

P(k) = k i k j g ij (k) 

where 

(&i(k)fi$(k)> 



S«( k ) = J™ ■ 

V — »oo 



V 



(6) 



(7) 



(8) 



and u(k) is the Fourier transform (FT) of the vector field 
u(q). Assuming that the latter is derived from a scalar 
potential as in Eq. (0) we have 



<7y(k) = kikjg(k) 



(9) 



where g(k) = Tr[g.y(k)] is a function of k = |k| only 
because the stochastic process is assumed to be statis- 
tically homogeneous and isotropic, and k = k/|k|. We 
thus have 



P(k) = k 2 g(k) = k 4 P${k) 



(10) 



where P<s>(k) is the PS of the fluctuations in the scalar 
field, i.e., 



P$(fc) = lim 

V->oc 



V 



(11) 



If one considers now a displacement field which varies 
as a function of time as in Eq. (0) , it follows that the PS 
of density fluctuations is proportional to the square of the 
function f(t). For a self-gravitating fluid such a behavior 
applies, and thus one can determine the function f(t) for 
the determination of the velocities 4 . Indeed Zeldovich 



3 We do not make the distinction here between physical and co- 
moving coordinates, and do not write the associated time depen- 
dent factors. Since we will analyze only IC for density fluctua- 
tions (and not velocities) these are not relevant details, and so 
we omit them for simplicity. 



4 Normally f(t) is chosen so that density perturbations are in the 
pure growing mode in which the velocity field is parallel to the 
displacement field. 



4 



originally proposed his approximation as an ansatz, on 
the basis that Eq. (|l|) implies the correct evolution of the 
density fluctuation in linearized Eulerian theory. The 
power of the ZA is that it can be applied well beyond 
the regime of validity of Eulerian perturbation theory, to 
which it matches at early times. 

To set up IC for the N particles of a cosmological NBS 
the procedure is then |32| : 

• to set-up a "pre-initial" configuration of the N par- 
ticles. This configuration should represent the fluid 
of uniform density pg. The usual choice is a simple 
lattice, but a commonly used alternative is an 
initial "glassy" configuration obtained by evolving 
the system with negative gravity (i.e. a Coulomb 
force) with an appropriate damping. 

• given an input theoretical PS Pth(k), the corre- 
sponding displacement field in the ZA is applied 
to the "pre-initial" point distribution. The cosmo- 
logical IC are usually taken to be Gaussian, and 
the displacements are determined by generating a 
realization of the gravitational potential 



*(q) = 2J flk cos (k ' l) + sin(k • q) 



with 



6 k = B, 



k 2 



(12) 



(13) 



where R\ and R2 are Gaussian random numbers of 
mean zero and dispersion unity. From Eq. (10) we 
see that this corresponds to generating a realization 
of a stochastic displacement field with PS <?ij(k) as 
in Eq. (0) and 



g(k) = P th (k)/k 2 



(14) 



III. 



ANALYTIC RESULTS IN fc-SPACE 



The configuration (or ensemble of configurations) gen- 
erated by the method outlined in the previous section has 
PS given through Eq. (|l^), and thus equal to the theo- 
retical PS Pth(k), up to the following approximations: 

• The system is considered as a continuous fluid. 
Thus we expect the exact PS of the (discrete) parti- 
cle distribution to differ by terms which come from 
the "granularity" (i.e. particle-like) nature of this 
distribution. 

• The calculations are performed at leading order in 
the amplitude of the density fluctuations, or equiv- 
alcntly, at leading order in the gradient of the dis- 
placements (cf. Eq. (||)). We would thus antici- 
pate that the exact PS of the generated configura- 
tions will have corrections which are significant for 
k larger than the inverse of a scale characterising 
the amplitude of the input PS. 



The rest of the paper is principally focussed on the 
consideration of the differences arising from the first 
point between the theoretical PS Pthik) and the exact 
PS (which we will simply denote P(k)) of the distribu- 
tion generated by the algorithm described in the previous 
section 5 . We refer the reader to |33[ |34| for analyses of 
the second point, i.e. of corrections coming from the use 
of the leading order ZA. These latter studies work in the 
continuum limit, and so completely decouple the problem 
of non-linear corrections from the effects of discreteness 
studied here. 



A. General results 

The starting point for our analysis is a result derived 
in jlj]. One considers, in d dimensions, the application of 
a displacement field u(r) to a generic point distribution. 
The latter is taken to have PS Pj n (k) and correlation 
function ^ n (r), given by the inverse Fourier transform 



&n(r) 



(27T) 



d d k, 



-ik-r 



(15) 



where the integral is over all space. The displacement 
field u(r) is assumed to be a realization of a continuous 
stochastic process, which is statistically homogeneous. 
An exact calculation Jl]] gives that the PS of the dis- 
tribution obtained in this way may be written as 



P(k) = J d d re-^p(k;r) (l +| m (r)) - (27r) d ,5(k). 



where 



p(k; r) 



d d i 



-ik-s 



p(s; 



(16) 



(17) 



and p(s; r) is the probability that two particles with a 
separation r undergo a relative displacement s. 

We note that our choice of notation here follows also 
that of |35| (rather than that of Q). In this work (see 
also Q) an expression for the PS generated by displace- 
ments given by the ZA is derived, for the case of a con- 
tinuous fluid. The general expression given is exactly 
that obtained by setting fi„(r) = in (|16|). This extra 
term in our expression arises because we do not make 
the approximation of treating the "pre-initial" configu- 
ration as a continuous uniform background. We note 
that this additional term contains not just the effect of 
taking into account the correlations in the "pre-initial" 
configuration, but also includes more generally all the ef- 
fects of the discreteness of the ("pre-initial" and final) 
distribution. In this respect we note that the correlation 



Note that the full PS is assumed to be a function of k, as it will 
not in general share the statistical isotropy and homogeneity of 
the theoretical PS (which makes it a function only of k = |k|). 
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function £in(r) for the "pre-initial" distribution contains 
generically a delta-function at r = 0, which is character- 
istic of its discreteness, as well as a non-singular function 
which describes correlations (for a detailed discussion see 
@@). 

From the definition of p(s, r) it follows that 

p(k, r) = J DuP({u(r)}) e -*M»M-»(<>)] ( 18 ) 

where the functional integral is over all possible displace- 
ment fields weighted by their probability P({u(r)}). For 
a displacement field which is (i) Gaussian, and (ii) sta- 
tistically isotropic (as well as homogeneous) it is then 
simple to show that 

p(k,r) = e - fafc ^«W (19) 

where a sum is implied over the labels i and j, and 

dy(r)=0y(O)-0«(r), (20) 

where 

9ij (r) = (ui(0) % (r)) (21) 

is the (matrix) two-displacement correlation function. 
We note that the scalar function 



fl (r) =Tr(ffy (r))=<u(0)-u(r)) 



(22) 



is simply the inverse FT of <?(fc) defined above (and, by 
statistical isotropy, a function only of r = |r|), and that 
g(Q) = (u 2 ) is the variance of the displacement field. We 
have that 

dij(r) = \([ui(0) - Ui(r)} [ Uj (0) - Uj (r)]) , (23) 

i.e., it is proportional to the correlation matrix of the 
relative displacements. 

Substituting (19) in ([l6|), we obtain 



P(k) = 



'krg-fcifc^-M (i + £ m ( r )) 



(27r) d «5(k). 



(24) 



It will be useful for our discussion to break this expression 
into two pieces, P(k) = P c (k) + Pd(k), written in the 
form 

P c (k) = J d d re-^ r (e^M^M - l) (25a) 
P d (k) = P in (k) + f d d re-^ ( e -*<M«(r) _ | m(r) . 

(25b) 



The first term Eq. ( 25a ) is the "continuous" piece of the 
generated PS (identical, as discus sed above to that given 
in |35)), and the second term Eq. (25b) is the contribution 
coming from the discreteness. 



B. Application to cosmological IC 

In the algorithm used to generate cosmological NBS, 
we have seen that the FT of gij{v) is [cf. Eqs. @ and 
@] given by 



9ij (k) 



Ph(fc) 
fc 2 



(26) 



Expanding the exponential factor in Eq. ( |25a| ) and ( p5 q ) 
in power series, we can thus obtain expressions for P c (k) 
and Pi(k) at each order in powers of P t h(k). At zero 
order we have evidently 



Pj°>(k)=0 
and, at linear order, 

P«(fc) = P th (fc) 



Pf (k)=P in (k) 



(27) 



(28a) 



k 2 



(2tt) 



d d q(k ■ q) : 



■Pth(q) 



X [P in (k + q)-P m (k)] (28b) 



To this order the PS of the generated distribution is thus 
the sum of the input theoretical PS and two discrete- 
ness terms: the PS of the "pre-initial" (i.e. lattice or 
glass) distribution and a second term which is a convo- 
lution of the input PS and the "pre-initial" PS. At next 
order in the expansion (i.e. at second order in Pth(fc)) 
we will obtain both further discreteness corrections, and 
corrections which survive in the limit in which we neglect 
discreteness completely. This result is in line with what 
we anticipated at the beginning of this section. 



C. Domain of validity of the expansion 

We have implicitly assumed above that the expansion 
we have performed is well defined 6 . This assumption cor- 
responds to that of finiteness of various integrals of the 
input PS Pth(fc). If the latter function is well-behaved, 
this corresponds to constraints on its asymptotic proper- 
ties, at small and large fc. To determine these constraints 
we consider a PS of the form 



P th (k) = Ak n f(k/k c ) 



(29) 



where A and n are constants, and f(x) is a function which 
interpolates between unity for a; < 1 and zero for x 3> 1, 
i.e., which may act as a cut-off for k > k c . In the use 
of this algorithm in cosmological simulations, for reasons 



6 We 
Eq. 
Eqs. 




that we have also assumed Gaussianity in deriving 
This is not in fact a necessary condition to obtain 
y the assumption that d;j(r) 
to recover the same result 



28b). Making instead o 
is bounded, it is easy (see also 
directly from an expansion of Eq. (I 
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which we will discuss further below, a very abrupt (usu- 
ally top-hat) such cut-off is always imposed at wavenum- 
bers of order the inverse of the scale characteristic of the 
interparticle separation 7 . Thus we will consider only the 
constraints at small k, i.e., the lower bound placed on the 
index n. 

Firstly we note that, using Eqs. (20) and (p6[), it is 
simple to show that (%(r) is well defined only if 



kmk d P th (k) =0, 



(30) 



i.e., if n > —d in Eq. (p9[). This is a condition which is 
always satisfied in cosmological models, as it follows from 
the finiteness of the one-point variance of the theoretical 
density fluctuations 8 . 

In App. ^ we analyse in detail the full expansion of 
P c {k) to all orders in A, separating two different cases: 
(i) — d < n < —d + 2, in which the one point variance 
g(0) is infinite, and (ii) n > — d + 2, in which case <?(0) 
is finite. From the expansions in each case one can infer 
the following: 

• For —d < n < —d + 2 the leading non-zero term, 
equal to Pth{k), approximates well the full P c (k) 
for the range of k in which 9 



A t \(fc) 



k d P th {k) « 1. 



(31) 



For — d+2 < n < 4 the criterion to satisfy the same 
condition is: 



at all orders 10 . In cosmological NBS the "pre-initial" 
distribution, as we have discussed, is usually taken to 
be a simple lattice or glass. We will see below that for 
the case of the lattice the coefficients of the expansion 
are sums which are regulated at small k, by the Nyquist 
frequency of the lattice (defined below). For the case 
of the glass, or indeed any distribution with an analytic 
Pj n (k), we limit ourselves to an analysis of the integral 
coefficient of the leading term in Eq. ( |28b| ). It is sim- 
ple to see, by Taylor expanding the expression inside the 
square brackets at small q, that the finiteness requires 
only the integrability of Pth(?) at small q. This coin- 
cides precisely with the condition Eq. (30). We expect 



Pj^k) will thus also be satis- 



that P d (k) » fg(k) 

fied when Eq. (|3l]) applies. We will verify below with 
numerical simulations that this is indeed the case. 

We note that the condition Eq. ( |3l| ) for the validity 
of the perturbative expansion at a given k is one which 
could be guessed from the simple continuum derivation 
using Eq. (g|), in which the expansion parameter is the 
amplitude of the theoretical density fluctuation: A 2 h (/c) 
is just a dimensionless measure of the amplitude of the 
density fluctuations in the theoretical model arising from 
wavenumbers around k 11 . Further we show in App. ^ 
that if condition Eq. (|}2|) is fulfilled for any k < k c , then 
Eq. ( |3l| ) is also. The two conditions are in fact essen- 
tially equivalent in the case that a cut-off is imposed as 
typically is done the cosmology. 



k*g(0) = k 2 {u z ) < 1 



(32) 



D. The leading non-trivial discreteness correction 



• For n > (—d+2) + d/2 there is, at next order in the 
expansion of P c (k), a correction proportional to k A . 
This implies that, for n > 4 the leading term Pth(k) 
is never well approximated at asymptotically small 
k by the PS of the generated IC. 

To analyse the expansion of the discreteness contribu- 
tion Pd(k) we need to specify the "pre-initial" distribu- 
tion. It is evident however that generically i t is a t least as 
convergent as than that of P c {k), since Eq. ( |25b| ) contains 
in the integrand simply an extra factor of (r) , which 
is typically smaller than unity and decreasing at large 
separations. For a Poisson distribution of number den- 
sity no, for example, one has £. n (r) = ^-S(r) (where S(r) 
is a Dirac delta function in <i-dimensions) , and therefore 
the expansion becomes trivial with Pd(k) = Pj n (fc) = 



' In this case the cut-off imposed in simulations, as explained be- 
low, is actually a function of k . 

8 The one point variance of density fluctuations is equal to £ t h(0), 
which is proportional [cf. Eq. ( |15[ )] to the integral of Pth(fc). 

9 In the cosmological literature Sr(k) is canonically defined with 
a numerical prefactor so that f(0) = (<5p 2 (0)} = / A 2 (fc)d(ln k). 
Given that the resultant factor depends on the dimension d we 
will not include it here. 



Let us now analyse in more detail the leading con- 
tribution to the generated PS arising from discreteness, 
i.e., the expression which we have denoted above by 
Pj 1} (k) 12 . We note [cf. Eqs. ( ]28a| - |28b| )] that this term 
arises at the same order as the input PS in the pertur- 



10 At leading order in the amplitude of the input theoretical PS P t ^ 
one therefore has P(k) = + P t f l (k). Thus for an exponent 

n < in (||) one will have P(k) m P th {k) for all k < (An ) 1/n . 
For n > 0, on the other hand, one can have P(k) fa Pth(k) 
at most in an intermediate range of k: at small k the Poisson 
variance of the "pre-initial" distribution will always dominate. 

11 Consistent with with Eq. this condition for the validity of the 
expansion can be stated equivalently in terms of the boundedness 
of the dimensionless quantity KKW-^WlK-W-^W])! > . ^ 
of the "gradient" of the displacement fields. We note that in a 
first version of the paper, a stronger condition was given for the 
validity of the expansion, n > — d + 2. This corresponds to the 
condition that variance of the displacement field be finite. While 
this stronger condition is assumed in the derivations in , and 
notably in arriving at Eq. (|19[), it is not a necessary condition 
for the validity of the method. We thank an anonymous referee 
for pointing out this error. 



In Appendix |b] 



we. 



.present some further analysis of the full ex- 
pansion of Eq. ( 25a ) , providing analytical expressions for some 
specific cases. 



l( L , , , ' 
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FIG. 1: Integral i(k) of Eq. (p7|), in three dimensions and av- 
eraged over shells of |k|, for an input PS Pth(k) as in Eq. ( psj ) 
with n — —2 and (i) f(k) — 1, i.e., without a cut-off (dotted 
line), and (ii) f(k) = 0(fczv — k), i.e., an abrupt (Heaviside 
step function) cut-off at k = fcjv (solid line). In the former 
case we see that I(k) is approximately constant for k < fcjv, 

and therefore the leading discreteness term Pj (fe) ~ k 2 for 
fc < fcjv- In the latter case, the discreteness term contributes 
only for k > fcjv. Such a cut-off is usually employed in the use 
of this algorithm in cosmological simulations. 



FIG. 2: As in previous figure, but now for n = and a top-hat 
cut-off function, for three different values of the cut-off. We 
see that as the cut-off increases the amplitude of I(k) does 
so (corresponding to the UV divergence of the sum Eq. j3^ ) 
for n > —1 in three dimensions). We see again that for a 
cut-off at fcjv the leading order discreteness term contributes 
only for k > fcjv, while for larger cut-off we have aliasing 
effects which manifest themselves in the appearance of the 
term P d (1) (fe) ~ k 2 for k < k N . 



bative expansion, i.e., at linear order in the amplitude of 
the input theoretical PS. We consider the specific case of 
a "pre-initial" distribution which is a simple cubic lattice. 
Its PS is 



P m (k) = (2^) d ^5(k-h) 

h^O 



where the sum over h is over all the vectors of the re- 
ciprocal lattice, i.e., h = m(27r/£), where I is the lattice 
spacing and m is a vector of non-zero integers. The min- 
imal value of |h| = 2tt/£, is the sampling frequency k s of 
the lattice, equal to twice the Nyquist frequency, which 
we will denote fcjv (and fcjv — 7r /l). It is instructive to 
rewrite the first order term Eq. (|28b|) in the form 



P^>(k) = Ak 2 k^- 2 I(li) 



(34) 



where 



n— 2 



x[P m (k + q)-P m (k)] (35) 

is dimensionless. Since Pm(k) = for k ^ h, we there- 
fore have, at linear order in our expansion in powers of 
the input PS, that, 

P(k) = P th (fc)+Pj 1} (k) 

= Ak n f{k/k c ) + Ak 2 k™~ 2 I(k) (36) 



where 



J(k) 



[k (h 



k)] 2 



h-kl : 



Ih-k| 



<N 



Ih-kl 



(33) for k ^ h. 



(37) 



The second (discreteness) term in Eq. (36) includes 
explicitly what is known as "aliasing" : power in the in- 
put spectrum at large wavenumbers (i.e. above the sam- 
pling frequency) gives rise to power at small k. Indeed 
the amplitude at small k of the discreteness term is pro- 
portional to I(k — > 0) in Eq. (|35|), which is a sum de- 
pending strictly on the power in modes at wave-numbers 
greater than or equal to the sampling frequency k s = 2fc jv ■ 
Further if one cuts at the Nyquist frequency k^, i.e., 
f(k/k c ) = 0(fcjv — k), where is the Heaviside step 
function, it follows that 7(k) = for k < fcjv. In this 
case therefore we have, for k ^ h, that 

P(k) = P th (k)e(k N -k) + Pi 1) (k)6(fc - k N ), (38) 

i.e., to leading order in the input spectrum the full PS of 
the generated IC is exactly equal to this input spectrum 
below the Nyquist frequency, and given by the discrete- 
ness term Eq. ( p5| ) above the Nyquist frequency. It is easy 
to verify that an analogous result applies if the cut-off is 
imposed in the first Brillouin zone (FBZ), i.e., setting 
the PS to zero but for vectors with all three components 
£ [— A; at, fcjv]- In cosmological simulations a cut-off is usu- 
ally imposed in this way (see e.g. ||, ps|). 

In Fig. [I] is shown the numerically computed value of 
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I(k) as a function of fc, in three dimensions 13 , for a pure 
power-law PS with n = —2, (i) without a cut-off [i.e. 
with /(fc) — 1 in Eq. (|2^)] and, (ii) with an abrupt top- 
hat cut-off at fcjy, i.e., f(k) — 0(kjy — k). In Fig. || we 
show the same quantity but for n = and three top-hat 
cut-offs at fcjv, 2fcjv, 3fcjv. We see clearly illustrated the 
behaviors discussed above. Note that for n < —d + 2 
(n > —d + 2) the expression for /(k) converges (diverges) 
without a cut-off, which explains the choices for the cut- 
off functions in the two figures. If a sharp cut-off is not 
implemented at fcjv wc scc that, in all cases, /(k) is non- 
zero and approximately constant for k < kjy. There is 
thus an associated aliasing term which is, to a very good 
approximation, proportional to k 2 . 



E. Accuracy of generation algorithm in k space 

We can now draw clear conclusions about the accuracy 
with which the generation algorithm, applied on a simple 
lattice, produces a point distribution with a PS approxi- 
mating the input PS of the form assumed in Eq. (p9|): 

• For — d < n < 4, and / an abrupt cut-off at fcjv, we 
have P{k) — P t h(k) for k < k N , up to corrections 
which depend parametrically on the dimensionlcss 
quantity A 2 h (fc) — k d P t h(k). For k > fcjy we have 

P{k) = P^p{k), where the latter is a discreteness 
term given explicitly in Eq. (j35|). 

• If any input power is included above the Nyquist 
frequency fc^ of the lattice (or, more precisely, out- 
side the FBZ of the reciprocal lattice), it leads to 
the appearance of power in the IC at fc < fcj\r (i.e. 
inside the first FBZ). With power included above 
the sampling frequency (fc = 2/cjv) there is an alias- 
ing term proportional to fc 2 at small k. In this case 
therefore the range of n which may be accurately 
represented at small k is limited to — d < n < 2. 

• For n > 4 one always has P(k) oc fc 4 at sufficiently 
small k, and the PS of the point process produced 
by the algorithm therefore does not approximate 
the input theoretical spectrum. 

• For n < —d the algorithm is not well defined be- 
cause the correlation function of relative displace- 
ments dij (r) is undefined. This is true in the infi- 
nite volume limit. In practice one generates IC in 
a finite system, usually taken to be a cube (with 
periodic boundary conditions). This means that 
in practice the input PS is always cut at the cor- 
responding fundamental frequency of the box, so 



We show the average for all vectors k with modulus in a bin 
centered about k = |k|. 



that, even for n < —d, the algorithm can be ap- 
plied. The implication of our result is that one will 
find in this case that the PS obtained will depend 
on this box size, becoming badly defined in the in- 
finite volume limit. We will verify that this is the 
case in our numerical study below. 



F. Glass pre-initial conditions 

The above conclusions were derived assuming that the 
"pre-initial" distribution is a simple lattice. The alterna- 
tive starting point quite often used in cosmological NBS 
are "glassy" configurations, obtained by evolving grav- 
ity with a negative sign and a strong damping on the 
velocities |||, p9fl . Without the damping, this system 
is essentially just what is known as the "one component 
plasma" in statistical physics (for a review, see ^]). The 
small k behavior of the power spectrum is then expected 
to be Pin(k) ~ k 2 at small k 14 . With the damping term 
what is found is a PS with a behavior ~ k A at small scales 
Assuming this form for the spectrum 15 it is easy 
to follow through the analysis given above for this case. 
The only change is that the term /(k) is now non-zero 
for all k : because Pi n {k) is non-zero for all k it is not 
possible to ha ve ze ro overlap of its support with that of 
Pth(fc) in Eq. ( |28b| ). This is what permitted this term to 
be zero in the case of a lattice and a top-hat cut-off at 
the Nyquist frequency. Thus in the case of a glass there 
will generically be a correction oc fc 2 at fc below the wave- 
number characteristic of the inter-particle distance in the 
glass. Thus the range of power-law spectra which may 
be accurately represented by the generation algorithm in 
this case at small fcis— d < n < 2. The models simu- 
lated in the context of cosmological N-body simulations 
are always well inside this range. 

What is the source of these limits on the representation 
of PS with n > 2 (or n > 4 on the lattice)? We remark 
that the appearance of such terms 16 would appear to 
be related to a well-known argument used by Zeldovich 
J42| , in determining the limits imposed by causality on 
fluctuations (See [fl4| and jl5) for discussion of this result 
and further references.): any stochastic process which 
moves matter in a manner which is correlated only up to 
a finite scale generates terms proportional to fc 2 in the PS 



Here "small" means compared to the inverse of the Debye length 
characterising the screening. This statement is true only if one 
neglects the damping, and assumes the system is in the fluid 
phase 

We assume thus that Pi n (k) ~ fc 4 up to k of order the "Nyquist 
frequency" (i.e. the inverse of a characteristic interparticle dis- 
tance) followed by a flattening to the required asymptotic form 
Pin(k) = l/ n o fo r la rger k. 

We note that in pa ] , which studies an input "top-hat" PS with- 
out power at small k, the fc 4 term in the PS has actually been 
actually measured numerically in the IC. The authors give it the 
same physical explanation we now discuss. 
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at small k. The coefficient of the k 2 term vanishes, leav- 
ing a leading term proportional to k A 1 if the additional 
condition is satisfied that the center of mass of the mat- 
ter distribution is conserved (i.e. not displaced) locally 
. The condition on the support of the displacement field 
required to make the coefficient of the k 2 vanish should 
thus be equivalent to a condition of local center-of-mass 
conservation under the effect of the displacement field. 



IV. RESULTS IN REAL SPACE 

We now turn to the consideration of the real space 
properties of the distributions generated by the algo- 
rithm. In this section we use the k space results of the 
previous section to determine these properties approxi- 
mately, but analytically. In the next section we will use 
numerical simulations in one dimension to show in detail 
the validity of these results. 



where W R {k) is the Fourier transform of the window 
function for a sphere of radius R, normalized so that 
W R (0) = 1. 

It is simple to show (see e.g. Q, and p^, Q for a 
more detailed discussion) that, for a PS of the form (|2S|), 
the behavior of the integral in (^2|) depends strongly on 
the value of n: 

• for —d < n < 1 the integral for <J 2 (R) is dominated 
by modes k ~ 1/R and one has 

a 2 (R) ~ k d P(k)\ k ^ 1/R « (43) 

• for n > 1 the integral is dominated by modes k ~ 
fcjT 1 (i.e. by the ultra-violet cut-off) and one has 
always 

° 2 ( R ) « -5^+1- ( 44 ) 



A. Definitions and background 

The quantities we will study in real space are the re- 
duced 2-point correlation function £(r) and the variance 
of mass in spheres. In fact we will principally consider 
the latter for reasons which we will explain below. 

We recall that £(r), for a statistically homogeneous 
distribution, is defined by 



(p(r)p(0)=p3(l + £(r-r')), 



(39) 



For n = 1 one obtains the transition behavior, in which 
the integral depends logarithmically on the cut-off k c . 
This gives <r 2 (R) cx \nR/R d+1 . 

The behavior in Eq. (|44|) is thus actually a limiting 
behavior. It is in fact a special case of a much more 
general result (see |57| |3^] for a discussion and references 
to the mathematical demonstration of this result): the 
most rapid possible decay in any mass distribution of 
the unnormalized variance of the mass ((AM) 2 )y in a 
volume V is proportional to the surface of the volume. 



where (...) is the ensemble average. For a discrete 
distribution (i.e. the case we always consider here) 
(p(r)p(r'))dVi dVi is the a priori probability to find 2 
particles in the infinitesimal volumes dV\, dVi respec- 
tively around ri and r 2 . The correlation function £(r) 
measures therefore the deviation of this probability from 
that in a Poisson distribution (equal to (po} 2 dVi dV^). It 
is related to the PS as its Fourier transform. 

The normalized mass variance o~ 2 {R) in spheres of ra- 
dius R is defined as 



a 2 (R) = 



(M(R) 2 ) - (M(R)Y< 
(M(R)) 2 



(40) 



where M(R) is the mass in a sphere of radius R, centered 
at a randomly chosen point in space. It is given in terms 
of the correlation function by 

^) = i/ d d rj d d r 2 £{\v x -v 2 \) (41) 

v V K ) Jv(R) JV(R) 

(where V(R) is the volume of a sphere of radius R), and 
in terms of the PS by 



o-\R) 



1 



(27r) d 



d d kP{\L)\W R {kt 



(42) 



B. Perturbative results in real space 

Returning now to Eqs. ( |27| - |28b| ), and using Eq.(|42"|), 
we infer that, at linear order in the amplitude of the input 
PS, we have 

a\R) = a 2 n (R) + a 2 h (R)+a 2 (R) (45) 
|(r) = &n(r)+£tfc(r)+&(r) (46) 

for the normalized mass variance and correlation function 
of the IC. The 'in' and 'th' subscripts in each case have 
the obvious meanings, with 'd' indicating the term asso- 
ciated to the linear order discreteness correction p9 (k). 
We have assumed implicitly that the integrals pick up 
negligible contribution from the regions, at large k, where 
the linear approximation to the full PS is not good. This 
will typically translate into a lower bound on R and r for 
the validity of Eqs. @ and @. 

It is simple to understand from Eqs. ( |45| ) and ( f46| ) why 
the question of the representation of real space properties 
of the IC generated using the ZA is non-trivially differ- 
ent from that of k space properties. In k space we had 
analogous expressions to Eqs. (^5|) and ([46|), from which 
it followed that P(k) w Pth(k) to very good accuracy at 
small k. One necessary ingredient for this was that the 
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term Pi n (k) could be neglected at small k, as it is iden- 
tically zero outside the FBZ on a lattice and decreasing 
very rapidly to zero (oc fc 4 ) in a glass. In real space we 
do not have the same "localization" at large k of the 
intrinsic fluctuations associated with the pre-initial dis- 
tribution. Indeed we have noted above that there is a 
limiting behavior (oc l/i? d+1 ) to the decay with radius 
R of the mass variance, for any distribution 17 . The am- 
plitude of this leading term is fixed by the inter-particle 
distance £, with of n ~ (£/R) d+1 , while that of the two 
other terms Eqs. (Eq) is proportional to the amplitude A 
of the input spectrum. Likewise for the correlation func- 
tion the intrinsic term £j„ (r) is generically delocalised in 
space, and depends only on the particle density, while the 
other two terms are proportional to the amplitude of the 
input PS. At sufficiently low amplitude, both quantities 
will be dominated at any finite scale by those of the un- 
derlying pre-initial point distribution, and thus will not 
be approximated by their behaviors in the input model. 
This is a behavior which is qualitatively different to what 
we have seen in reciprocal space. We now examine in a 
little more detail these two-point quantities. We treat 
them separately as they are quite different for what con- 
cerns their comparison to the continuous theoretical in- 
put quantities: being an integrated quantity, the mass 
variance is intrinsically smooth and can be directly com- 
pared with its counterpart in the input model. 

C. Mass variance in spheres 

Given Eq. (|4q) , and the limits we have discussed on the 
behavior of the variance, we can immediately make a sim- 
ple classification of the PS of the form ( p9| ) for what con- 
cerns the representation of their variance in real space. 
The faithfulness of such a representation requires simply 

a* h (R) » *l(R). (47) 

For either a lattice or glass we have the "optimal" de- 
cay af n (R) oc l/R d+1 . In order for Eqs. (fi|) and @) 
to be valid we require that Eqs. ( p^ - |28b| ) be valid. As 
discussed in the previous section we expect this to corre- 
spond to the criterion that A 2 h (k) = k d P t h{k) < 1 for the 

relevant k. Given that P^ 1 (k) is at most proportional to 
k 2 at small k, the associated variance is also oc \/R d+1 
above the interparticle distance i, and thus sub-dominant 
with respect to the leading term at all scales. Since we 
generically cut the input spectrum around fc/v, and will 
consider simple power law spectra up to this scale with 
n < —d, it suffices to have 

A 2 N = A t 2 h (fc N ) - k d P th (k N ) < 1 . (48) 



Up to a numerical factor of order unity this is none other 
that the criterion 18 that <J 2 h (£) < 1, and so it follows 
that we expect the following behaviors: 

1. For n > 1 we have seen that cr 2 h (R) ~ 1/R d+1 , i.e., 
aj h (R) has the same functional behavior as that of 
the "pre-initial" variance. Given that the former 
is necessarily smaller at the inter-particle distance, 
the condition Eq. ( f47| ) will never be fulfilled, as 
the full variance will be dominated by that of the 
pre-initial configuration. 

2. For -d < n < 1 we have that a 2 h {R) ~ l/R d+n , 
which thus decays more slowly than the "pre- 
initial" term. Thus there will be a scale R m in such 
that for R > R m i n one can satisfy the condition 
Eq. (|47j). It is easy to infer that, for any d, we 
have 

2 

Rmin ~ I (-^A (49) 

D. Two point correlation function 

The case of the two point correlation function is simi- 
lar. The determination of the range of faithful represen- 
tation of the theoretical correlation function is, however, 
more complicated by the very non-monotonic behavior of 
the correlation function in both the (unperturbed) lattice 
and glass. This leads, as we will explain, to a strong de- 
pendence on how the correlation function is smoothed 
when it is estimated in a sample. 

Unlike for the variance, there is no intrinsic limit on the 
rapidity of the decay of the correlation function for point 
processes. Indeed for a Poisson process one has £in( r ) = 
for r > 0, and exponentially decaying correlation func- 
tions are commonplace in many physical systems. For 
both a lattice and glass distributions the leading term 
£i n (r) in Eq. ( f46| ) presents a very non-trivial behavior. 
The two point correlation function of the lattice is in fact 
not a function of r, but a distribution which depends on 
r: it is proportional to a Dirac delta function when r links 
any two lattice points, and equal to —1 otherwise (see Ap- 
pendix |^ for the explicit expression) . For the glass the 
correlation function is not known exactly — it depends 
on the details of the generation of the glass configura- 
tion used — but generically it will be expected to have 
a similar oscillatory structure describing its very ordered 
nature, with decay only at scales considerably above the 
interparticle distance 19 . This underlying highly ordered 



While the result we cited concerning the variance applies strictly 
to the case of statistically homogeneous and isotropic distribu- 
tions, it can be shown (see [m], f38|) that it applies also to the 
variance in spheres measured in a lattice. 



For the case n > 1, this is true only because the input PS is cut 
at the Nyquist frequency; for n < 1 it is true even without the 
cut-off. 

The characteristic property of these configurations is that the 
force on particles is extremely small. This imposes a very strong 
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structure is evidently not washed out by the application 
of very small displacements. In particular for relative 
displacements much smaller than the initial interparti- 
cle separation, it is clear that the form of the underlying 
correlation function will remain highly oscillatory up to a 
scale considerably larger than the interparticle distance. 
Just as in the case of the mass variance, therefore, one 
can conclude that the theoretical term in Eq. (|4^) will al- 
ways be dominated by the discreteness terms up to some 
scale, which becomes larger as the input amplitude is 
decreased. 

A simple analytical estimate, like that given above for 
the variance, of the scale at which the theoretical term 
will dominate the discreteness terms, and thus at which 
the input theoretical two-point correlation function is 
well approximated by that of the generated IC, is not 
possible: for the lattice such an estimate must take also 
into account the term £d(r) which together with £m(r) 
gives a regular oscillating and decaying function; for a 
glass we do not have the analytical form of the correla- 
tion function. 

There is a further important difficulty if one wishes to 
compare the correlation function in generated IC with 
the input one. In estimating the correlation function 
in a finite sample one must introduce a finite smooth- 
ing: one computes it by counting the number of pairs of 
points with separations in some finite interval, typically 
a radial shell of some chosen thickness. Indeed while 
the full correlation function is in general a function of r, 
this procedure makes it a function of r like the theoret- 
ical correlation function. Given that, at low amplitude 
of the relative displacements, £(r) has both a strongly 
oscillating and strongly orientation dependent behavior, 
such a smoothing can change very significantly its behav- 
ior. Thus the scale at which agreement may be observed 
between the measured ensemble averaged two point cor- 
relation function and that of the input model will depend 
both on A, I and the precise algorithm of estimation of 
the correlation function. 



NUMERICAL STUDY IN ONE DIMENSION 




FIG. 3: PS of a model with n = -1/2, sharp cut-off f(k) = 
Q(k - fcjv) and Ajv = 1.77 x 10~ 3 (A = 1(T 3 ). The simu- 
lation results are averaged over one thousand realisations of 
IC, generated using the standard algorithm (adapted to one 
dimension) . 



real-space mass variance for a large ensemble of config- 
urations, which is not numerically feasible (for modest 
computational power) in three dimensions. The exact 
ensemble average results given above, on the other hand, 
are easily calculated. The simplified and more explicit 
expressions for the relevant quantities are given in Ap- 
pendix |^. There is no intrinsic difference of importance 
between one and three dimensions for the questions we 
address 20 . 

We consider the case in which the pre-initial distribu- 
tion is a lattice. Following our discussion in the previous 
sections we study separately the four following specific 
examples for input PS as in Eq. (29): (i) n=-l/2 (exam- 
ple of — d < n < 1), (ii) n=3 (example of 1 < n < 4), 
(iii) n=7 (example of n > 4) and (iv) n=-2 (example of 
n < — d, in which case we have found the algorithm to be 
badly defined in the infinite volume limit). We will spec- 
ify the cut-off function in each case. We then also present 
numerical results for the two point correlation function in 
just the first of these models to illustrate the discussion 
of this quantity given at the end of the preceding section. 



In this section we study the generation algorithm us- 
ing numerical simulations. This allows us to verify our 
conclusions about two point properties in reciprocal and 
real space, derived in the limit of small amplitudes of the 
input PS. Further it allows it to show the accuracy of the 
full analytic expression Eq. (p4[), for any input amplitude. 
We work in one dimension because of the numerical fea- 
sibility of the study in this case: we measure directly the 



A. n = -1/2 (Case -d < n < 1) 

In Fig. U are shown results for an input PS Pth{k) = 
Ak^ 1 / 2 with A = 10~ 3 , which corresponds to Ajv = 
1.77 x 1CP 3 . Here, as in the rest of this section, we 
use units of length in which the inter-particle distance 
is equal to unity. We have imposed a sharp FBZ cut-off 
f(k) = Q(k — kjy). In the figure we see, as expected, ex- 



correlation between the positions of particles. In studies of the 
one component plasma, mentioned above the appearance 

of multiple peaks in the correlation function is observed as the 
temperature is lowered. 



One minor exception for the case of the two point correlation 
function, related to the last point discussed in the previous sec- 
tion, is discussed at the appropriate point below. 
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FIG. 4: PS of a model with n = —1/2, continuous cut-off 
f(k) = exp(-fc/2fc]v) and Ajv = 1.77 x 10~ 3 (A = 1CT 3 ). 



FIG. 6: Mass variance in spheres of a model with n = —1/2, 
sharp cut-off f(k) — 0(fc — fcjv) and two different amplitudes 
of the theoretical PS. 




FIG. 5: PS of a model with n = -1/2, sharp cut-off f(k) 
Q(k - fcjv) and Ajv = 0.18 (A = 0.1). 



cellcnt agreement between the PS measured by averaging 
over a thousand realisations of IC, generated using the 
standard algorithm (with Gaussian displacements) in a 
periodic interval containing a thousand particles, and the 
theoretical expression at linear order, as given in the pre- 
vious section. Note that on the x — axis is given /c/2/cjv, 
so that first Bragg peak appears at unity, and the sharp 
change in the PS at 0.5. 

Fig. ^ differs only in that we have now imposed a con- 
tinuous cut-off / = e - k / 2k N _ Again we observe, as ex- 
pected, excellent agreement between the measured PS 
and the theoretical expression. The agreement between 
the input PS and the measured PS is, however, less per- 
fect around fcjv, because the discreteness term P9~\k) 
contributes now inside the FBZ (i.e. for k < fcjy). The 
effect is, however, very small as the latter term is, in this 
range, proportional to k 2 . 

In Fig. H are shown results for the same shape PS, but 
now with a higher amplitude, A — 0.1, corresponding 
to Ajv = 0.18. The cut-off here is sharp. Shown are 
the input theoretical PS, the average over one thousand 



realisations, and the exact expression for the PS. We are 
not in this case in the regime in which the perturbative 
expansion of the full PS is valid at k^, and therefore do 
not plot P^\k) as in the previous figures. Indeed we see 
that the PS of the generated IC begin to deviate sensibly 
from the input theoretical IC already at a k significantly 
smaller than k^ , with a discrepancy of about a factor of 
two in the amplitude at k = /cat. Note that, nevertheless, 
the results of the simulations agree extremely well with 
the exact expressions for the full PS. 

In Fig. ^ are shown the real-space variance in spheres 
of radius R (i.e. intervals of length 2R) for the case of 
the sharp cut-off and the two different amplitudes just 
considered. The curves labelled "exact" are those cor- 
responding to the ensemble average of the full IC, and 
those labelled "theoretical" are those of the input model. 
We see clearly illustrated the results anticipated in the 
previous section: for low amplitudes the exact curve is 
dominated at small distances by the variance of the un- 
derlying lattice, and the low amplitude theoretical ex- 
pression (which has a behavior a 2 (R) <x 1/R 1 ^ 2 ) is ap- 
proximated only once this term coming from the lattice 
(with <J 2 (R) oc l/R 2 ) has decayed sufficiently. At the 
higher amplitude the theoretical expression, on the other 
hand, is well approximated for scales just above the in- 
terparticle distance 21 . 



B. n = 3 (Case 1 < n< 4) 

Figs. to [l(] show exactly the same quantities as the 
four previous figures, but now for an input power-law PS 



21 The discrepancy between the variance appears smaller than that 
in the PS (shown in Fig. |jj ) at the inverse scale due to the 
different range of scale on the y-axis in the two plots. The relative 
difference is in fact of the same order. 
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FIG. 7: PS of a model with n = 3, sharp cut-off f(k) 
Q(k - fcjv) and Ajv = 0.35 {A = 1.8 x 10" 5 ). 



FIG. 9: PS of a model with n = 3, sharp cut-off f(k) 
Q(k - k N ) and Ajv = 35 (A = 0.36). 




with n = 3. The two amplitudes chosen are given in the 
captions, the low amplitude corresponding to the case 
where the linear approximation to the exact formula for 
the PS is a good approximation. Figs. ^ to || illustrate 
the more important difference that arises in the case that 
n > 2 when the cut-off imposed on the PS is smooth in- 
stead of being imposed sharply inside the FBZ: the k 2 
term at small k generated in P^\k) in this case dom- 
inates the input PS at small k so that it is no longer 
faithfully represented by the PS of the generated IC at 
any k. Fig. ^ shows essentially the same thing as Fig. [|. 
For higher amplitudes the agreement of the input PS with 
that of the generated IC is shifted to smaller k. The ex- 
act formula for the PS agrees very well with that of the 
generated IC measured in the simulations, but the lin- 
ear approximation to the discreteness effects at larger k, 
given by P^\k), is no longer a good approximation. 

Comparison of Fig. |Io| with Fig. ^ shows the difference 
between the cases n < 1 and n > 1 for what concerns 
the behavior of the mass variance in real space. Because 
the theoretical variance has the same scale dependence 



(a 2 (R) cx 1/i? 2 ) as the lattice variance, the latter always 
dominates the former if the amplitude is low. Specifically, 
if the input mass variance at the lattice spacing is less 
than that of the lattice (which is of order unity) the mass 
variance of the IC is not approximated at any scale by 
that of the input model. 



n = 7 (Case n > 4) 



Figs. 11 



to [jj] show results for the PS of a single low 
amplitude n — 7 input model, for the case of a sharp and 
continuous cut-off respectively. These figures illustrate 
the limitation discussed in the previous section for the 
representation of a small k input PS with n > 4. Us- 
ing the sharp cut-off inside the FBZ the term P^' (k) is 
zero for k inside the FBZ, but nevertheless the theoret- 
ical behavior at small k is not represented because the 
corrections to Eq. ( |38| ) , at quadratic order in the ampli- 
tude A, are non-zero. Thus at asymptotically small k we 
see the PS of the generated IC deviate from the input 
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FIG. 11: PS of a model with n = 7, sharp cut-off f(k) 
6(fc - fcjv) and An = 1.76 x 1(T 3 (A = 1.85 x 1CT 5 ). 



FIG. 13: PS from simulations of a model with n = —2, sharp 



cut-off f(k) = 0(k— fcjv) and Ajv 
for different number of particles. 
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FIG. 12: PS of a model with n = 7, continuous cut-off f(k) — 
exp(-fc/2fciv) and Ajv = 1.76 x 10~ 3 (A = 4.8 x 1(T 5 ). 



one 22 . Further the behavior at the smallest k is well fit 
by a fc 4 behavior, which is shown in Appendix ^ to be 
that of the quadratic order correction. 

In Fig. |l2] one can observe the dramatic effect, as we 
saw illustrated also in Fig. ^, of using a continuous cut- 
off for the case n > 2. Just as in the case of n — 3 we 
see that the input PS is no longer well approximated — 
indeed not even poorly approximated — by the PS of the 
generated IC. Note that, differently from Fig. O, there is 
no range of intermediate k where the input PS is approx- 
imated. This is because it is the correction P^p (k) which 
dominates at small k, with an amplitude proportional to 
same (linear) power of the input PS. Correspondingly in 
Fig. [ll] the k at which a deviation towards the k 4 behav- 



22 The numerical integration of the exact expression in this case is 
very difficult because of a very rapidly oscillating behaviour in 
d(x) at large x. The 'exact' curve has thus been calculated just 
far enough at small k so that the deviation from the input PS 
may be discerned. 



ior is observed can be shifted to arbitrarily small k by 
taking a sufficiently low initial amplitude. 



D. 



n = -2 (n < -d) 



We show finally in Fig. [13] results for the PS for aver- 
ages over simulations of the case n = —2. In this case, 
as we have discussed above, the algorithm is not well de- 
fined in the infinite volume limit, because the variance of 
relative displacements at any scale is a divergent. The 
implementation of the algorithm in a finite sample, with 
periodic boundary conditions, is perfectly well defined as 
the spectrum of modes is cut-off at small k by the fun- 
damental, fixed by the box size. In the figure we show 
the results for the PS of the averages of 1000 generated 
configurations, for different numbers of particles, i.e., for 
different sizes of the system. As anticipated the results 
depend strongly on the box size, and neither the ampli- 
tude nor the shape of the input PS is approximated well 
by that of the generated distributions. 



E. Two point correlation function 

Figs. |l4| and [l^ illust rate q uantitatively the discussion 
and conclusions in Sec. 



[VP above. They show both the 



exact two point correlation function, and a smoothing of 
it, for IC corresponding to an input power-law PS with 
n = —1/2. The smoothing is defined by a convolution of 
the discrete density distribution with a spatial window 
function Wl- 



Pc(r) 



cPr'W L (\T-T'\)p d (T'), 



(50) 



where p c ( r ) is the density function of the continuous field, 
Pd(r) of the discrete distribution and L is the characteris- 
tic scale introduced by the smoothing. For the correlation 
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FIG. 14: The absolute value of the measured and theoretical 
two point correlation function for a theoretical input model 
with power-law PS and n = —1/2, for amplitude A — 2 x 



icr 



The curve labelled 'exact' is the result of a numerical 



evaluation of the full expression of the correlation function 
for the given model. The curve 'smoothed' gives the same 
quantity, but now smoothed by a convolution with a Gaussian 
window function as given in Eq. © with W L (k) = e" fc2 . 
The 'theoretical' curve is the correlation function of the input 
model, proportional to 1/r 1 ' 2 . 
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FIG. 15: The same quantities as in Fig. |l4| for a much larger 
amplitude, A — 0.2 of the input PS. The smoothing is exactly 
the same as in the previous figure. 



function this gives 



dr'FT^, [\W L (k)\ 2 ]i(x'), 



(51) 



where FT denotes the inverse Fourier transform of 
Wi(fc). For the latter we have taken here a simple Gaus- 
sian form as specified in the caption of the figures. 

We observe that the range of agreement between these 
quantities and the theoretical correlation function is dif- 
ferent — illustrating that the result depends on the 
smoothing — and, further, that this range depends also 
on the amplitude of the input model. Just as for the 
mass variance, the scale above which the theoretical and 



measured quantities converge increases (for any given 
method of estimation/smoothing) as the amplitude de- 
creases. Further for sufficiently low amplitude pertur- 
bations the underlying structure of the lattice becomes 
visible if the estimated two point function is resolved to 
the required level (by a sufficiently narrow smoothing). 

One remark is appropriate here on the relation between 
these results and those in three dimensions. One impor- 
tant effect in that case is not illustrated by these results: 
if one takes, as is usually done, a simple pair estimator 
for £ (r) using spherical shells of equal width, the volume 
of the shells grows as r 2 . Therefore the oscillations of the 
true lattice or glass correlation function will be attenu- 
ated much more rapidly as a function of distance than by 
the smoothing considered here in one dimension. This, 
however, does not change any of the conclusions above: 
the scale at which agreement will be observed between 
the measured and theoretical quantities will depend on 
the size of the bins, and taking sufficiently small bins one 
can always make the oscillatory structure of the under- 
lying correlation function dominate for a sufficiently low 
amplitude of the input spectrum. 



VI. SUMMARY AND CONCLUSIONS 

We first summarize our findings on the accuracy and 
limitations of the standard algorithm for generating IC 
for cosmological simulations. We then discuss the con- 
clusions we can draw, in the light of our analysis, about 
the some numerical work on IC |2(| which partly mo- 
tivated our study. Finally we turn to the relevance of 
our results to the problem of understanding discreteness 
effects in the evolution of cosmological simulations. 



A. Results on generation algorithm 

We have investigated systematically the algorithm 
used to generate IC of N-body simulations in cosmol- 
ogy, for any given input PS. More specifically we have 
focussed on the comparison of the two point correlation 
properties, in real and reciprocal space, of the IC with 
those of the input theoretical models. We consider input 
PS which are a simple power-law P(k) oc k n , but the 
corresponding results for more complicated cases may be 
easily inferred. Our main results are: 

1. Applied on a grid with appropriate sharp cut-off at 
the Nyquist frequency fc^v, the point distribution 
produced by the algorithm has PS exactly equal 
to the input one, below k^, to linear order in the 
amplitude and for — d < n < 4. For k > k^ we 
have also given the exact expression for the PS, 
which is thus the leading discreteness correction in 
this space. It is a term of high amplitude, with a 
damped oscillating form with maxima at the Bragg 
peaks of the underlying lattice. 
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2. Applied to a 'glass' pre-initial configuration, the 
result is almost the same, except that the discrete- 
ness correction has a small k tail proportional to k 2 . 
Thus the range of "faithful representation" of the 
PS is— d < n < 2. This latter restriction is not of 
relevance to current cosmological models, for which 
the effective exponent at all k is within this range. 

3. The algorithm does not produce IC representing 
faithfully an input PS with n > 4 for arbitrarily 
small k. There is the case because there is a term 
proportional to k 4 in the PS of the generated PS, 
at second order in the amplitude of the input PS. 

4. For the case n < —d the algorithm is not well de- 
fined in the infinite volume limit, and we have veri- 
fied that results in a finite volume depend strongly 
on the volume. 

5. The transposition of these results to real space is 
more subtle than one might have anticipated, due 
to the fact that the mass variance and two point 
correlation of the underlying 'pre-initial' point dis- 
tribution are delocalised in this space. 

6. For models with —d<n< 1 the real space variance 
in spheres can be well represented by the generated 
configurations starting from a finite scale R m in pro- 
portional to the inter-particle spacing. For typical 
chosen input amplitudes it is a few times this dis- 
tance, but we note that it diverges as the amplitude 
of the input model goes to zero. 

7. For models with n > 1, the real space variance is 
always dominated, at linear order in the amplitude, 
by the "pre-initial" variance of the lattice or glass. 

8. The conclusions concerning the representation of 
the reduced two point correlation function arc quite 
similar to those for the mass variance: the theo- 
retical properties are recovered above a finite scale 
proportional to the inter-particle distance, which 
diverges as the amplitude goes to zero. In practice 
there is a further difference with respect to the mass 
variance, in that the value of this scale depends 
also on the smoothing is necessarily introduced in 
estimation of the correlation function. For a suffi- 
ciently narrow smoothing the correlation function 
will always show at a given scale, for sufficiently 
low amplitude of the input model, the underlying 
structure of the lattice or glass configuration. 



B. Comments on precedent literature 

Let us now consider, in the light of these results, the 
articles |2l], |2^, |2j| which have partly motivated this 
work. These two collaborations draw, on the basis of 
numerical studies, very different conclusions about the 



measured mass variance in spheres and two-point corre- 
lation function of the IC of cosmological NBS. 

In cosmology the IC of NBS are invariably studied only 
in reciprocal space, simply because it is the natural one 
for the description of cosmological models at early times. 
In the first of these papers |2C]] the authors examined 
instead IC in real space, through a numerical study of 
the IC of some large cosmological simulations performed 
by the Virgo consortium Their finding was, very 

surprisingly, that the measured and theoretical values of 
both the mass variance in spheres and the two point cor- 
relation function did not match. In Q| the same analysis 
was repeated by a different set of authors, and an error 
in the normalization in ppf of the theoretical variance 
was identified. Correcting for this error the authors con- 
cluded that the agreement between the measured and 
theoretical properties was good for the variance, while 
the authors of [^(| , in a reply [^2| , argued that the agree- 
ment was still very poor. For the two point correlation 
function the results of both sets of authors agreed, show- 
ing an estimated correlation function qualitatively and 
quantitatively different to the expected one. The two 
sets of authors gave a quite different interpretation to 
this discrepancy: in [ pp[ it was attributed to a probable 
systematic difference between the two quantities due to 
the underlying correlation in the "pre-initial" configura- 
tion, while pl| argued that it was more likely simply due 
to statistical noise of the estimator. In a further article 
p3| the second authors analyzed these same quantities 
in the IC of another set of cosmological simulations, and 
arrive at the same conclusions as in pi) concerning both 
quantities. 

For what concerns the mass variance we have seen 
that the degree of agreement between the theoretical and 
measured variance depends on the normalization of the 
model, i.e., on the initial red-shift of the simulation. Nei- 
ther collaboration has studied the dependence of their 
conclusions on this crucial parameter, nor identified it 
as relevant. Thus the conclusions of |2l], |2^] about the 
reliability in general of the representation of the input 
mass variance by the IC are, strictly, incorrect. However, 
their conclusion that the representation of this quantity 
is good for the specific set of IC considered — normalized 
at an amplitude which is typical in practice in cosmolog- 
ical simulations - is correct. That is the agreement they 
observe in a modest range (see e.g. the figure 3 in p3[). 
from a few times the interparticle distance to a scale ap- 
proaching the box size, at which finite size effects start 
to play a role, is real (rather than purely accidental as is 
implicitly suggested by |2(], |22| ) . However the dominant 
lattice term can clearly be identified at smaller scales, 
and it is evident in view of our discussion that the range 
of agreement will decrease (and ultimately disappear) if 
one considers the same model with a lower normalization. 

For the two point correlation function, we have seen 
that the degree of agreement depends not only on the 
amplitude of the input model, but also on the details of 
the spatial smoothing in the estimator. Again neither 
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collaboration has pinpointed explicitly the importance of 
this consideration in evaluating the faithfulness of the 
representation. The authors of |2(], ^2|), however, are 
correct when they argue that the difference observed is 
a systematic one, and that the oscillating behavior ob- 
served is due to the correlations in the underlying pre- 
initial (lattice or glass) configuration. In attributing the 
difference to "noise" the other group is incorrect, insofar 
as such noise would be a finite size effect which should 
disappear in the ensemble average. However, noise can 
of course play a crucial role in a finite sample in hiding 
the underlying signal in the regime in which it may, in 
principle, approximate well the theoretical model. 



C. Physical relevance of results on IC 

We have considered in this paper solely the question 
of the accuracy with which the standard algorithm for 
generating IC for cosmological NBS represents the theo- 
retical correlation properties. This question is essentially 
interesting insofar as it is relevant to the question ad- 
dressed by the series of articles of which this is the first: 
the quantification of the differences between the results 
of evolved N body simulations and the corresponding the- 
oretical quantities. This question will be addressed fully 
in the subsequent papers, and we limit ourselves now to 
a partial discussion of the physical relevance of our find- 
ings. 

The most important result from a practical point of 
view is that, at linear order in the theoretical density 
perturbations, there is a contribution to the PS of the IC 
additional to the theoretical PS. This is a source for grav- 
itational structure formation through the Poisson equa- 
tion, which in a given simulation cannot be separated 
from the theoretical term. Indeed we note that the linear- 
ity of this term in the amplitude of the relative displace- 
ments means that, if the early time evolution follows the 
Zeldovich approximation, this term is amplified linearly, 
just like the theoretical term. On the other hand, it con- 
tributes significantly only above the Nyquist frequency, 
and therefore, given that gravity tends to transfer power 
very efficiently from large to small scales (see, e.g., @)j 
one would expect its effects to be washed out over time. 
However if one wishes to quantify precisely discreteness 
effects, our quantification of this leading discreteness con- 
tribution in the IC is an important first step. 

In quantifying such effects it is important also to first 
understand the recovery of the continuum limit. Our 
results here, as we will now discuss, actually are quite 
informative in this respect. Let us consider the limit 
in which one recovers exactly the properties of the the- 
oretical continuum model. Given an input theoretical 
model for a cosmological NBS, we introduce two param- 
eters with the standard method of discretization we have 



discussed here 23 : £, the lattice spacing in physical units, 
and the initial red-shift Zi (which fixes the amplitude A 
of the input PS, with A — > as z\ — > oo). 

The continuum limit should evidently correspond to 
taking I — > (and thus kjq — ► oo). Let us consider first 
taking I — ► at fixed z\. This corresponds in our anal- 
ysis above to working at fixed amplitude of the PS. Our 
results above tell us that the representation of the PS 
is good provided we satisfy the condition Eq. (^Tj) for 
the validity of the perturbative expansion. This quan- 
tity in fact converges to zero for k^ 3> k c , and so the 
criterion for good agreement in k space for all k is sim- 
ply A 2 h (/c c ) -C 1. This agreement becomes arbitrarily 
good as we take Z\ — > oo (i.e. z\ — > ). Likewise in real 
space, it follows from Eq. (|4^) that we converge towards 
an arbitrarily good representation of the mass variance 
when the limit is taken in this way. The same is true of 
the two point correlation function. Thus the correlation 
properties of the discretized IC converge exactly to the 
continuum IC. 

Our results concerning the differences in real space 
quantities concern the limit Z\ — > oo, at fixed I. We 
have seen that there is, in this case, no convergence to- 
wards the continuum model. Thus, in the IC, the order 
of the limits in i and z\ cannot be interchanged. It will be 
shown in the companion paper p4j , that the same non- 
commutativity of the limits is observed in the evolved 
systems. This in fact is just a specific example of a well- 
known fact about the validity of continuum Vlasov dy- 
namics to describe a system with long range interactions 
|^|, |^] . In this context it is known and well documented in 
certain systems that the continuum limit is approached 
as N — » oo keeping the time of evolution fixed, while 
taking the time to infinity first one diverges from the 
collisionless limit (see, e.g., pq]). Lowering the initial 
amplitude of a NBS increases the time of evolution (up 
to a given time), and thus the behavior we are inferring 
from the analysis of the IC corresponds to this same one. 

These comments on the continuum limit are also of 
practical relevance, as they tell us how one should study 
convergence to this limit numerically (in order to un- 
derstand the precision of results). It follows from what 
we have just discussed that it is best to keep z% fixed 
as the particle density is increased. Further the contin- 
uum limit can only be defined clearly in the presence 
of a cut-off in the input PS, with the continuum limit 
being approached when the interparticle distance is de- 
creased well below the inverse of this scale. In most of 
the numerical studies in the literature on discreteness in 
cosmological NBS these points have not been taken into 
account 24 . Indeed we note that the very widely used 



In reality there is of course also the box size L, which we have 
taken in our study to be infinite. The finite particle number N 
is given by (L/£) 3 . 

An exception is some of the cited work of Melott et al.. Some sets 
of simulations are compared in which only the particle density is 
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standard software package COSMICS for generating IC 
p9| fixes automatically the initial red-shift of the simu- 
lation when the physical particle density is given, rather 
than leaving it as a free parameter, making such con- 
trolled tests difficult. Indeed if no cut-off is imposed in 
the input PS, the criterion used to fix the red-shift makes 
it increase with the particle density. These points will be 
further discussed in forthcoming work. 

We recall finally that our results on the limitations of 
the use of the algorithm for very blue spectra are of rel- 
evance to some studies in the literature of gravitational 
evolution from such spectra. Specifically we note their 
usefulness in understanding quantitatively results in J2{| 
and |2(|. These studies consider gravitational N body 
simulations (in two and three dimensions, respectively) 
starting from IC generated on a lattice using the stan- 
dard algorithm discussed here, taking input theoretical 
PS with vanishing initial power in some range of small 
k: in [^H) a top-hat PS is used, while in [£6| a gaus- 
sian centred on a chosen wavenumber. In both cases our 
results show that there is a term proportional to k in- 
duced at small k already in the IC, which will dominate 
at small k. The explicit expression for this term, which 
arises at second order in the expansion of t he co ntinuum 
piece P c (k) of the full PS, is given in Eq. (|AlJ). In @ 
the dominant contribution from the k A term in the IC at 
small k is observed numerical ly, an d indeed the authors 
relate it (as discussed in Sect. [II F above) to Zeldovich's 
argument about "minimal power". In p6 |, on the other 
hand, the k 4 term is seen (and observed, as expected, 
to grow with an amplification proportional to the square 
of the linear growth factor) only after some time. The 
authors describe in this case the fc 4 tail as "generated" 
by the dynamical evolution, which is evidently not quite 
accurate as the term is in fact present (albeit at lower 
amplitude) already in the IC. 

We are indebted to Andrea Gabrielli for extensive 
discussions and explanations of his results reported in 
We thank Thierry Baertschiger and Francesco Sy- 
los Labini for numerous useful conversations, and Alvaro 
Dominguez and Alexander Knebe for helpful comments 
on the first version of this paper. We are indebted to an 
anonymous referee for pointing out an important error in 
the first version of this paper. 



APPENDIX A: PROPERTIES OF THE 
EXPANSION OF P c (k) 

In this appendix we study in more detail the pertur- 
bative expansion used in the paper of P c (k), Eq. (2- r );i). 

To simplify our analysis we will take the function dy (r) 
[defined in Eq. (|20|)] to be diagonal and isotropic, i.e., 



varied, keeping both the initial amplitude and the cut-off in the 
input PS fixed in units of the box size. 



dij{v) = d(r)5ij. This allows us to obtain simple analyt- 
ical results, which are exact in one dimension and which 
we expect to be valid only with minor modifications in 
three dimensions. 



Expanding Eq. (25a) in powers of d(r) we have 



P c (k) = V(-fc 2 )" 1 / d d re- lk - r [d(r)Y 



(Al) 



We will suppose a theoretical PS in the form of Eq. (|2£ 
(and assume that g'y(k) = Sijg(k) = 5ijPth{k)/k 2 ) 



1. Case -d < n < -d + 2 

We can work in this case without the UV cut-off in 
the PS, since d(r) is well defined without it [cf. Eqs. p0[ . 
Their evaluation gives 

d(r) = — r(n-l)sin(^r 1 " n [d=l](A2) 
= lr(n) S in(^Vi-" [d = 3] (A3) 



The integrals in (Al) are then divergent as r — > oo, but 
defined in the sense of distributions. Evaluating them we 
obtain 



P c (k) 
where 



m— 1 



p(™)(fc) = 



rn—1 



A r - 



ik rn(n+d)-d _ (A4) 



-A- 



2tt 



- — sin —rmrfn — 1) r(l + m — mn) 
ml \ 2 



(r(n-l)sin(^))" 1 [d = l] 



.4: 



-r(2-m(l + n)) 



-m(n + l)?r j (r(n) sin (^j 



(A5) 
(A6) 
[d = 3] 



We note that the expansion (A4) is in fact an asymp 



totic expansion, i.e., it is strictly divergent, but if an 
appropriate finite number of terms are taken, for any 
given fc, it approximates closely the well defined function 
Eq. (25a). This behavior is sho wn i n Fig. [l6|, in which 
is plotted the ratio of the series (A4) summed up to the 
771-th term, and P c (k). We see that for the ratio first 
converges to unity as m increases, but then diverges at 
progressively smaller k for m > 30. 

P t h{k) is well approximated by P c (k) if 



Ak n > A 2 a 2 k 



1{n+d)-d 



i.e., for 



Aa-ik 



n+d 



< 1. 



(A7) 
(A8) 
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FIG. 16: Ratio of the series (A4) summed up to the m-th 
term to the exact (numerically evaluated) P c (k). We observe 
that there is first clear convergence (i.e. approach to unity) 
as m increases, and then divergence at larger values of m. 



It can be checked using Eqs. flAq) or (A6) tha t a m /a m _i 
is of order unity for small m, so that Eq. (A8) corre- 
sponds to the criterion Eq. (pll) given in the paper. Note 
that we can rewrite Eq. ( |A8]) in terms of the variance of 
mass in spheres of the theoretical fluctuations, using the 
approximation (e.g. (3J 



38|) 



a 2 (R) = bk d P th (k)\ k= 



R- 1 - 



(A9) 



where the coefficient b is of order unity. The condition 
for faithful representation of the PS of the input model 
at wavenumber k can thus be written: 



a 2 (R)\ 



fc=K-i < I- 



The case — d + 2 < n < 00 



(A10) 



In this case we must include the UV cut-off in the PS, 
in order that d(r) be well defined. The latter is then not 
a simple power-law at all scales as in the precedent case, 
and we ar e un able to compute analytically the terms of 
the series (Al). We can, however, compute very simply 
the first corrections to P t h{k). Since g(0) is finite (for 
—d + 2 < n < 00), we can rewrite Eq. (Al) as 



(All) 



where we have used the identity 
{2ir) d S(k) = (2Tr) d e- k29{0) 5(k) = e- k2g(0} [ d d re 



ik-r 



(A12) 

Exp andin g first the exponential containing g(r) in 
Eq. (All) we obtain 



P c (k) = e- k s(0) 
x (k 2 g{k) + k 4 



d d r[g(r)] 



2 e -,k. r 



(A13) 
0(k e [g(r)f) 



Expansion of the exponential pre-factor then gives 

P c (k) ~ P th (k) + (A14) 
k 4 [ I d d r[g{r)Ye-^ r - g(0)g(k) 



By d imen sional analysis one can see that the integral in 
Eq. ( A14) scales as (c\ + C2k d+2n ~ A ) , where c\ and C2 are 
non-zero constants. For the range of index n considered 
it follows that: 

• For — d + 2 < n < 2, the dominant correction to 
Pth comes from the term oc g(0)g(k) and therefore: 



P c (k) ~ k*g{k) - k*g(0)~g(k)- 



(A15) 



It follows that the condition for a faithful repre- 
sentation of the theoretical PS (Pth(k) = k 2 g(k)) 
is 



Mk 2 « 1, 



(Ai6) 



which corresponds to the condition Eq. ([Il"|). For a 
sharply cut-off theoretical PS 



Pth{k) = 



Ak n for k < k c . 
otherwise. 



one has 



5(0) 



Akl 



tt(ti — 1) 
2n 2 {n + 1) 



[d=l] 
[d = 3] 



(A17) 

(A18) 
(A19) 



Dropping the nume rical factors, for simplicity, the 
condition Eq. ( Al6| ) can be written as 



n+l 



[d=l] 

[d = 3] 



(A20) 

(A21) 
(A22) 



Since we are consideri ng th e case n > —d + 2 here, 
this means that Eq. (A16) is, for k < k c , a more 
restrictive criterion than that found in the previ- 
ous case. However, since Aj? h (£;) is a monoton- 
ically increasing function of k up to k c , the two 
conditions are essentially equivalent in cosmologi- 
cal NBS, in which one generically imposes a cut-off 
around fc^y. We note further that the condition is 
then also equivalent to 



a 2 (R)\ kN=R -i < 1, 



(A23) 



which is equivalent to (A10) (since cr 2 (R) is in this 
case also a monotonically decreasing function of R) . 
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• For 2 < n < 4 the main correction comes from the where 



integral in Eq. ( |A14| ) 



Pc(k) 



k 2 ~g(k) + -fc 4 



d d r[g(r)] 2 e 



(A24) 



For the sharply cut-off theoretical PS of Eq. flAlTp , 
the integral can be evaluated analytically in the 
limit k — > 0. This gives 



P c (k) 



Ak n 
Ak n 



A' 



,fc 4 k 



2n-3 



A 



2tt 2n - 3 
fc 4 fc 2 "- 1 
2^2 2n- 1 



[d = 1] (A25) 
[d = 3] . (A26) 



Up to numerical factors of order unity the leading 
correction is the same as in the previous case, and 
thus the same criteria apply for the validity of the 
perturbative expansion as in the previous case. 

• For n > 4 the resulting PS is dominated by the fc 4 
correction. The full expression for P c (k) therefore 
does not approximate P t h{k) at sufficiently small 
k. 



k N n, 



(B5) 



and n are triple integers. The smallest q in the sum (B4) 
is the Nyquist frequency, so that th e lead ing term at small 
k is the one we discussed in Sect. Ill D| , proportional to 
fc 2 . 



In one dimension all the terms in the series (B4) may 
be calculated analytically, for the case — d < n < —d + 2 
without a UV cut-off. The leading fc 2 term is 



Pd(k) 



2Afc"- 2 C(2 - 



n)fc 2 +C(fc 3 ). 



(B6) 



We can then estimate the scale fc up to which this term 
is sub-dominant compared to the input spectrum: 



k<, (2C(2-n)) 



l/n-2 



(B7) 



We note that this scale is independent of the amplitude 
of the PS (and for n = -1/2, fc<4.2fcjv). This result 
is completely in line with the numerical calculations of 
the contributions of these terms for various other cases 



presented in Sect. ([IIP) 



APPENDIX B: 



DISCRETENESS CORRECTIONS 
TO THE PS 



We analyse further in this appendix the full expansion 
to al l ord ers of the exact discreteness correction in the PS 
Eq. ( [25b| ). Then for the specific case of an input PS with 
—d < n < —d + 2, and no UV cut-off, we can evaluate 
the expression analytically in one dimension. This gives, 
in particular, an analytic expression for the coefficient of 
the leading contribution, proportional to fc 2 and allows 
a precise determination of the range of fc in which this 
term is sub-dominant with respect to the input PS. 

Expanding the exponential in Eq. (25b) we have 



AP d (k) = P d (k)-P in (k) (Bl) 

CO „ 

= E(- fc2 r/ d d re-^ r [d(x)} m ^ n (r)- 



m—1 



which can be rewritten as 



APd(k) = (2^F £ ( " fc)2m L A# (m) (qm»(q + k) 

(B2) 



where 



L>< m >(k) := FT[(d(x) m ], 



where FT denotes the FT as defined in Eq. 
pre-initial simple cubic lattice this gives 

oo 

AP d (k)= ^(-fc) 2m ^£>( m )(q + k) 
m=l q^O 



(B3) 
For a 

(B4) 



APPENDIX C: ANALYTICAL RESULTS IN ONE 
DIMENSION 

In this appendix we give some simplified analytic ex- 
pressions for the exact PS, mass variance and two point 
correlation function in one dimension. We have made use 
of these expressions in our numerical study of various in- 
put PS in Sect. 0. 

We recall first the correlation properties of a simple 
cubic lattice (in d dimensions for generality) which we 
will take as the "pre-initial" distribution in what follows. 
For the reduced two point correlation function one has 

6«t(ri,r 2 ) - (p(r)p(r')) - 1 = ^8 (n - r 2 - 1) - 1, 

( C1 ) 

where 1 is a generic displacement vector of the lattice. 
The expression Eq. ( |33"| ) is simply the Fourier transform 
of this expression. 

Let us now consider the case of one dimension. To 
compute the variance we use its expression as a function 
of the PS (see p§): 



— 1 ( + °°dk<- 



Z7T 



P(k) (C2) 



kR J 

or, equivalcntly, as a function of the correlation function: 



1 



' (R) " 8* 



dx^(x) x 



(C3) 



x [-2x0(x) + (x- 2R)9{x - 2R) + (x + 2R)9{x + 2R)} , 



where 6(x) is the Heaviside function. Using Eqs. (|C2j) 
or (P3[) with (|33"|) or (CI) respectively, we obtain the 
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following result for the variance of a lattice with grid 
spacing equal to unity : 



+ 00 

E 



sin(27rmi?) 
2nmR 



(C4) 



As anticipated in the previous section we obtain the same 
limiting behavior of the variance at large scales as for a 
homogeneous and isotropic distribution with PS P(k) ~ 
k n and n > 1, i.e., a 2 (R) ~ 1/R d+1 with d=l. 

We now compute an expression for the PS directly from 
( p4| ) , for the case of a one-dimensional system and a "pre- 
initial" lattice configuration. Using Eq. (CI) and rear- 
ranging terms we obtain: 



+ OC 



P(k) = exp(-fc 2 ff (0)) ^(fc - 2ttZ) 



(C5) 



1= — 00 



where d(x) = q( 0) — The first term on the right 



hand side of Eq. (C5) contains all the divergent terms in 
the PS. The second term is a regular function of k which 
has the behavior P(k) ~ k 2 g(k) at small k if g(k) ~ fc a 
with a < and P(fc) ~ fc 2 if a > 0, unless J^t^oo S (0 = 
0, in which case P(k) ~ k 2 g(k) also for a > 0. 

Performing a Fourier transform of Eq. we obtain 
the correlation function in the form 



1 

2^ 



+00 



d{x') 

x (i + e m K)) -i- 



g-^-x') 2 /^^') x 



(C6) 



Note that in the limit that no displacements are applied 
(i.e. d(x) — > 0), the argument of the integral is S (x — x'). 
Thus we recover explicitly for small displacements £ (x) ~ 



£ m (x) + . . . . Substituting Eq. (|Clj) in Eq. (|Cj) we then 
obtain the result for the specific case of a "pre-initial" 
lattice configuration: 



+00 



£» = -!+ J2 



l=—oo 



4wd(l) 



(C7) 



To obtain the vari ance we use the same procedure. Using, 
for example, Eq. (|C^) with Eq. we get: 



a 2 (i?) = -1 + 



dx (1 + £in(x))V djx) X 



x [h{x,2R)+h(x,-2R)-2h(x,0)] 



+ 00 



— / dx(l+£ in (x)) x 



(C8) 



J2 e- lkl [exp(~~k 2 d(l)))-exp(-k 2 g(0))}, where 



x [-2f(x,x) + f(x-2R,x) + f(x + 2R,x)] 



f(x, y)=x erf \^—j=j , /t(x, y) = e . (C9) 

Expanding at small d(x) it is possible to obtain also ex- 
plicitly an expression of the form <J 2 (R) — v 2 at (R) + .... 
In the specific case of an initial lattice distribution the 
variance can be written: 



i= — 00 

x [h(l,2R)+h(l,-2R)-2h(l,0)] 



(CIO) 



^ 2 [-2/(1, i) + /(i - 2J2, + /(i + 2i2, i)] • 



8R 2 
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